O’Hara on GitHub 

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')

source(here('common_fxns.R'))

1 Summary

Create a species richness map for threatened marine species.

Combine taxa-level trend maps to create maps of number of threatened species overlapping intensifying or deintensifying stressors, per cell.

Plot both by absolute numbers of species impacted in a location and by proportion of species (impact counts / species richness).

This is all done at the level of cumulative impact category, not individual stressors.

2 Methods

2.1 Create species richness map for this assessment

Set up the dataframe of species to include. By taxon, load all species range maps, smash down to taxon-level species richness; then combine these for a total species richness.

spp_incl <- get_incl_spp()

spp_sens <- spp_incl %>%
  filter(!is.na(stressor))
cell_id_rast <- raster(here('_spatial/cell_id_mol.tif'))
ocean_a_rast <- raster(here('_spatial/ocean_area_mol.tif'))

spp_rich_rast <- raster(here('_output/rasters/n_spp_map.tif'))
values(spp_rich_rast)[values(spp_rich_rast) == 0] <- NA

2.2 Setup for mapping across all taxa

For each impact category, sum increasing and decreasing stressor locations for count:

  • gather all taxa-level map files for increasing or decreasing stressors into a raster stack, by level of intensification
  • collapse the stack using calc(fun = sum, na.rm = TRUE)

2.2.2 Plot net intensification

Net intensification is the difference between the number of species under intensification and the number under abatement.

str_cats <- c('fishing', 'ocean', 'land-based', 'climate', 'all')

infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s%s.tif')
  ### stressor category, incr/decr

for(str_cat in str_cats) {
  # str_cat <- str_cats[5]
  for(level in 1:3) {
    incr <- raster(sprintf(infile_stem, str_cat, 'incr', level))
    decr <- raster(sprintf(infile_stem, str_cat, 'decr', level))
    
    diff <- calc(stack(incr, -decr), fun = sum, na.rm = TRUE)
    
    diff_pct <- diff / spp_rich_rast
    diff_pct[incr == 0 & decr == 0] <- NA
        
    hcl.pals('diverging')
    brks <- seq(-1, 1, .05)
    clrs <- hcl.colors(n = length(brks), palette = 'Red-Green', rev = TRUE)
    plot(diff_pct,
         col = clrs,
         breaks = brks,
         axes = FALSE,
         main = sprintf('Net intens, Pct spp, %s strs level %s', str_cat, level))
  
  }
}

2.2.3 Plot intens vs abate

Because a zero net intensification can result from low intensification cooccuring with low abatement, or high of each, let’s examine the frequency with which various combinations occur. Just examine level 1 as this is the most likely to have overlaps.

str_cats <- c('fishing', 'ocean', 'land-based', 'climate', 'all')

infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s1.tif')
  ### stressor category, incr/decr

for(str_cat in str_cats) {
  # str_cat <- str_cats[5]

  incr <- raster(sprintf(infile_stem, str_cat, 'incr'))
  decr <- raster(sprintf(infile_stem, str_cat, 'decr'))
  
  df <- data.frame(incr = values(incr), 
                   decr = values(decr), 
                   nspp = values(spp_rich_rast)) %>%
    filter(!is.na(nspp)) %>%
    filter(!is.na(incr) & incr != 0 & !is.na(decr) & decr != 0) ### we're interested in co-occurrence
  ### only 57989 cells where increase and decrease co-occur?
  
  df1 <- sample_frac(df, size = 1)

  tmp <- ggplot(df1, aes(x = decr, y = incr)) +
    geom_bin2d(bins = 15) +
    geom_abline(slope = 1, intercept = 0, color = 'red') +
    scale_fill_viridis_c() +
    scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50, 100, 200)) +
    scale_y_log10(breaks = c(1, 2, 5, 10, 20, 50, 100, 200)) +
    labs(title = sprintf('Intens vs abating, %s strs level 1', str_cat),
         x = 'nspp abating', y = 'nspp intensifying')

  print(tmp)
}

# install.packages("remotes")
# remotes::install_github("slu-openGIS/biscale")

library(biscale)
library(cowplot)

ocean_a_df <- data.frame(prop = values(ocean_a_rast),
                         cell_id = 1:ncell(ocean_a_rast)) %>%
  mutate(a = prop * (res(ocean_a_rast)[1])^2 / 1e6)

### skip ocean - only one stressor, no spp can simultaneously be increasing
### and decreasing
str_cats <- c('fishing', 'land-based', 'climate', 'all')

infile_stem <- here('_output/rasters/intens_maps/intens_%s_%s.tif')
outfile_stem <- here('figs/intens_biplot_%s.png')

for(str_cat in str_cats) {
  # str_cat <- str_cats[4] ### 4 = all
  outfile <- sprintf(outfile_stem, str_cat)
  
  if(!file.exists(outfile)) {
    message('Creating bivariate plot for ', str_cat, ': ', basename(outfile))

    incr <- raster(sprintf(infile_stem, str_cat, 'incr'))
    decr <- raster(sprintf(infile_stem, str_cat, 'decr'))
    
    df <- data.frame(incr = values(incr), 
                     decr = values(decr), 
                     nspp = values(spp_rich_rast),
                     cell_id = 1:ncell(incr)) %>%
      filter(nspp > 0 & !is.na(nspp)) %>%
      mutate(incr_pct = incr / nspp,
             decr_pct = decr / nspp)
             
    ### create classes: three of intensification (incr), and three of abatement (decr).
    ### types: "quantile" (default), "equal", "fisher", and "jenks".
    ### apparently use "fisher" instead of "jenks" for larger datasets
    data <- bi_class(df, x = incr_pct, y = decr_pct,
                     style = "fisher", dim = 3) %>%
      mutate(bi_class = factor(bi_class),
             bi_int = as.integer(bi_class))
    ### using jenks:
    #     1-1     1-2     1-3     2-1     2-2     2-3     3-1     3-2     3-3 
    # 2068375   47836   25942 1058722   14297    1768  257262    2743      99 
    ### using fisher:
    #     1-1     1-2     1-3     2-1     2-2     2-3     3-1     3-2     3-3 
    # 2140894   39477   26031 1005903   12893    1254  248652    1841      99 
    ### Quantiles fails b/c non-unique values (in the function).
    # quantile(df$incr_pct, seq(1/6, 1, 1/6))
    # quantile(df$decr_pct, seq(.05, 1, .05)) ### little of the ocean is decreasing!
    # incr_cutoff <- quantile(df$incr_pct[df$incr_pct != 0], .25)
    # decr_cutoff <- quantile(df$decr_pct[df$decr_pct != 0], .25)
    # data <- df %>%
    #   mutate(x = case_when(incr_pct == 0 ~ '1',
    #                        incr_pct <= incr_cutoff ~ '2',
    #                        incr_pct <= 1 ~ '3',
    #                        TRUE ~ 'error'),
    #          y = case_when(decr_pct == 0 ~ '1',
    #                        decr_pct <= decr_cutoff ~ '2',
    #                        decr_pct <= 1 ~ '3',
    #                        TRUE ~ 'error')) %>%
    #   mutate(bi_class = paste(x, y, sep = '-'),
    #          bi_class = factor(bi_class),
    #          bi_int = as.integer(bi_class))
    ### table(data$bi_class)
  
    bi_pcts <- data %>%
      left_join(ocean_a_df, by = 'cell_id') %>%
      mutate(a_tot = sum(a)) %>%
      group_by(bi_class) %>%
      summarize(a_level = sum(a),
                pct_level = round(a_level / first(a_tot), 4)) %>%
      separate(bi_class, into = c("x", "y"), sep = "-") %>%
      mutate(x = as.integer(x), y = as.integer(y))
      
  
    
    bi_class_levels <- levels(data$bi_class)
  
    bi_map <- map_to_rast(data, 'bi_int') %>%
      rasterToPoints() %>%
      as.data.frame() %>%
      mutate(bi_level = bi_class_levels[layer])
    
    # create map
    intens_map <- ggplot(data = bi_map, aes(x, y, fill = bi_level)) +
      geom_raster(show.legend = FALSE) +
      bi_scale_fill(pal = "DkViolet", dim = 3) +
      ggtheme_map() +
      labs(title = sprintf('Intensification and abatement: %s', str_cat))
  
    intens_legend <- bi_legend(pal = "DkViolet",
                        dim = 3,
                        xlab = "+ % intens",
                        ylab = "+ % abate",
                        size = 8) +
      geom_text(data = bi_pcts, aes(x, y, label = pct_level),
                size = 3, color = 'white') +
      scale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) +
      theme(panel.background = element_blank(),
            plot.background = element_blank())
      
    # combine map with legend
    final_plot <- ggdraw() +
      draw_plot(intens_map, 0, 0, .8, 1) +
      draw_plot(intens_legend, 0.6, .6, 0.4, 0.4)
  
    ggsave(plot = final_plot, filename = outfile, height = 4, width = 6, dpi = 150)
  }
  knitr::include_graphics(outfile)
}